NMDS on understory vegetation in Coastal Barrens group (MA & LI)
Kathleen Stutzman
2024-02-26
Libraries:
# libraries -------------------
library(ggplot2)
library(vegan)
library(ggvegan)
library(tidyverse)
library(ggordiplots)Set working directory
knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')Functions:
# functions -------------------Run NMDS
# 3 dimensions
mali_3d <- metaMDS(rel_mali5,
distance = 'bray',
k = 3,
trymax = 20,
autotransform = FALSE)## Run 0 stress 0.1272962
## Run 1 stress 0.1275739
## ... Procrustes: rmse 0.03533478 max resid 0.1292549
## Run 2 stress 0.1272962
## ... New best solution
## ... Procrustes: rmse 0.0009770007 max resid 0.003222686
## ... Similar to previous best
## Run 3 stress 0.1275739
## ... Procrustes: rmse 0.03479693 max resid 0.1269719
## Run 4 stress 0.1275743
## ... Procrustes: rmse 0.03490449 max resid 0.1272927
## Run 5 stress 0.1272965
## ... Procrustes: rmse 0.001136387 max resid 0.0037354
## ... Similar to previous best
## Run 6 stress 0.1275741
## ... Procrustes: rmse 0.03482774 max resid 0.1271008
## Run 7 stress 0.1272958
## ... New best solution
## ... Procrustes: rmse 0.0003339852 max resid 0.00103124
## ... Similar to previous best
## Run 8 stress 0.1272962
## ... Procrustes: rmse 0.0006479442 max resid 0.002175514
## ... Similar to previous best
## Run 9 stress 0.1272961
## ... Procrustes: rmse 0.0005806935 max resid 0.001964567
## ... Similar to previous best
## Run 10 stress 0.1272962
## ... Procrustes: rmse 0.0006353023 max resid 0.002144382
## ... Similar to previous best
## Run 11 stress 0.1272958
## ... Procrustes: rmse 0.000133377 max resid 0.0003272971
## ... Similar to previous best
## Run 12 stress 0.127296
## ... Procrustes: rmse 0.0005204546 max resid 0.001740775
## ... Similar to previous best
## Run 13 stress 0.127574
## ... Procrustes: rmse 0.03500575 max resid 0.1278142
## Run 14 stress 0.1275736
## ... Procrustes: rmse 0.03486528 max resid 0.1273871
## Run 15 stress 0.1272961
## ... Procrustes: rmse 0.000635704 max resid 0.002127893
## ... Similar to previous best
## Run 16 stress 0.1275742
## ... Procrustes: rmse 0.03505738 max resid 0.1279639
## Run 17 stress 0.1275737
## ... Procrustes: rmse 0.03491215 max resid 0.1275222
## Run 18 stress 0.1275735
## ... Procrustes: rmse 0.0348122 max resid 0.1272177
## Run 19 stress 0.1275739
## ... Procrustes: rmse 0.03497939 max resid 0.1277474
## Run 20 stress 0.127296
## ... Procrustes: rmse 0.0002188625 max resid 0.0005793351
## ... Similar to previous best
## *** Best solution repeated 8 times# converges, stress of 0.127 - 0.128
# 2 dimensions
# mali_2d <- metaMDS(rel_mali5,
# distance = 'bray',
# k = 2,
# try = 40,
# trymax = 40,
# autotransform = FALSE)
# stress still too high to represent in 2 dimensions, lowest is ~ 0.2Significant species & envr factors for Axis 1 & 2
## NMDS1 NMDS2 Species pval
## GAPR -0.43283432 -0.36592469 GAPR 0.026
## VAAN -0.58548007 0.23659024 VAAN 0.004
## QUCO 0.47191797 0.31542190 QUCO 0.010
## PIRI -0.09464334 0.63995157 PIRI 0.001
## KAAN -0.64346682 0.04951010 KAAN 0.001
## MELI -0.21626012 -0.54646216 MELI 0.003
## RUSP -0.62458049 0.06725700 RUSP 0.004
## SMGL 0.06322859 0.74451448 SMGL 0.001
## QUPR -0.58237657 -0.08846107 QUPR 0.011## NMDS1 NMDS2 env.variables pval
## Time_from_Treat 0.3664852 -0.3459599 Time_from_Treat 0.046
## Moss -0.1485437 0.5470263 Moss 0.013
## BA_HA 0.5192040 -0.3902518 BA_HA 0.006
## PIRI.BA_HA 0.5939817 -0.4423917 PIRI.BA_HA 0.002Significant species are: GAPR, VAAN, QUCO, PIRI, KAAN, MELI, RUSP, SMGL, QUPR
Significant environmental factors are: Moss, Basal area per hectare, and PIRI basal area per hectare.
I’m only going to keep one BA measure; PIRI has a lower p value, so thinking that one, but open to other ideas.
3D Plot
This is a 3D ordination being represented by 2D
Using the first two axes, which represent the most variation in the ordination space
GG plot of Axis 1 vs. Axis 2
Plot with significant species & envr factors Axes 1 & 2
Plot with indcator species & envr factors: Axes 1 & 2
Plotting envr factors & species for axes 1 & 3
## NMDS1 NMDS3 Species pval
## QUIL 0.001238461 0.6996833 QUIL 0.002
## GAPR -0.413939366 0.3911541 GAPR 0.019
## VAAN -0.571941725 -0.2857426 VAAN 0.003
## QUCO 0.436621926 -0.4459518 QUCO 0.005
## KAAN -0.615704312 -0.2950916 KAAN 0.001
## RUSP -0.572890031 -0.4210261 RUSP 0.001
## QUAL 0.163594982 -0.5134256 QUAL 0.020
## QUPR -0.561433797 -0.2528658 QUPR 0.004## NMDS1 NMDS3 env.variables pval
## avgLD 0.2099211 0.6656505 avgLD 0.003
## BA_HA 0.5201249 0.3122623 BA_HA 0.014
## PIRI.BA_HA 0.6201426 0.2041855 PIRI.BA_HA 0.003Significant species are: QUIL, GAPR, VAAN, QUCO, KAAN, RUSP, QUAL, QUPR
Significant environmental factors are: average litter depth, Basal area per hectare, and PIRI basal area per hectare.
Graph of Axis 1 vs. Axis 3: indicator species
Graph of Axis 1 vs. Axis 3: significant species
Plotting envr factors & species for Axes 2 & 3
## NMDS2 NMDS3 Species pval
## QUIL 0.03601922 0.6998788 QUIL 0.001
## PIRI 0.62635710 -0.2909393 PIRI 0.001
## MELI -0.52040825 -0.2251256 MELI 0.008
## QUAL -0.32982027 -0.5095111 QUAL 0.003
## SMGL 0.74047454 0.1515597 SMGL 0.001## NMDS2 NMDS3 env.variables pval
## Moss 0.5336551 -0.1538819 Moss 0.018
## Mineral_Soil 0.4287884 -0.3477778 Mineral_Soil 0.013
## Wood 0.1803424 0.4688395 Wood 0.043
## Litter_Duff -0.2796757 -0.4292763 Litter_Duff 0.041
## avgLD -0.1887942 0.6465422 avgLD 0.002Significant species are: QUIL, PIRI, MELI, QUAL, SMGL
Significant environmental factors are: moss, mineral soil, wood, litter-duff, & average litter depth
Graph of Axis 2 vs Axis 3: indicator species
Graph of Axis 2 vs Axis 3: significant species
Run PerMANOVA to look at the impacts of region and treatment on community composition
summary(as.factor(ma.meta.data_veg$Treat_Type))## Control FallRx Harvest MowRx SpringRx
## 5 3 6 3 6ma.tt <- adonis2(rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)
print(ma.tt)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 4 1.9569 0.31903 2.1082 0.003 **
## Residual 18 4.1769 0.68097
## Total 22 6.1339 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pairwise follow up
library(pairwiseAdonis)
pairwise.adonis2(rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)## $parent_call
## [1] "rel_mali5 ~ Treat_Type , strata = Null , permutations 999"
##
## $FallRx_vs_SpringRx
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.21549 0.10956 0.8613 0.596
## Residual 7 1.75138 0.89044
## Total 8 1.96687 1.00000
##
## $FallRx_vs_MowRx
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.48859 0.46035 3.4122 0.1
## Residual 4 0.57276 0.53965
## Total 5 1.06134 1.00000
##
## $FallRx_vs_Harvest
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.40379 0.20243 1.7767 0.121
## Residual 7 1.59094 0.79757
## Total 8 1.99473 1.00000
##
## $FallRx_vs_Control
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.39033 0.23712 1.8649 0.107
## Residual 6 1.25581 0.76288
## Total 7 1.64614 1.00000
##
## $SpringRx_vs_MowRx
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.65697 0.28336 2.7678 0.011 *
## Residual 7 1.66151 0.71664
## Total 8 2.31848 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $SpringRx_vs_Harvest
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.4952 0.15598 1.8481 0.042 *
## Residual 10 2.6797 0.84402
## Total 11 3.1749 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $SpringRx_vs_Control
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.16729 0.0666 0.6422 0.772
## Residual 9 2.34457 0.9334
## Total 10 2.51186 1.0000
##
## $MowRx_vs_Harvest
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.72529 0.32578 3.3823 0.017 *
## Residual 7 1.50107 0.67422
## Total 8 2.22636 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $MowRx_vs_Control
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.84147 0.41918 4.3303 0.019 *
## Residual 6 1.16594 0.58082
## Total 7 2.00741 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Harvest_vs_Control
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 1 0.58418 0.21102 2.4072 0.038 *
## Residual 9 2.18412 0.78898
## Total 10 2.76830 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## attr(,"class")
## [1] "pwadstrata" "list"Significant: SpringRx vs. MowRx; SpringRx vs. Harvest (borderline at 0.051); MowRx vs. Harvest; MowRx vs Control; Harvest vs Control
Indicator species analysis
library(labdsv)
ma.meta.data_veg$Treat_Group <- NA
ma.meta.data_veg$Treat_Group <- ifelse( ma.meta.data_veg$Treat_Type == "FallRx", 1, ma.meta.data_veg$Treat_Group)
ma.meta.data_veg$Treat_Group <- ifelse( ma.meta.data_veg$Treat_Type == "MowRx", 2, ma.meta.data_veg$Treat_Group)
ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "SpringRx", 3, ma.meta.data_veg$Treat_Group)
ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "Harvest", 4, ma.meta.data_veg$Treat_Group)
ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "Control", 5, ma.meta.data_veg$Treat_Group)
mali_isa <- indval(mali5_avg, ma.meta.data_veg$Treat_Group)
summary(mali_isa)## cluster indicator_value probability
## QUPR 2 0.8549 0.003
## KAAN 2 0.7534 0.010
## RUSP 2 0.7097 0.008
## SMGL 4 0.6896 0.012
##
## Sum of probabilities = 4.667
##
## Sum of Indicator Values = 6.68
##
## Sum of Significant Indicator Values = 3.01
##
## Number of Significant Indicators = 4
##
## Significant Indicator Distribution
##
## 2 4
## 3 1gr <- mali_isa$maxcls[mali_isa$pval <= 0.05]
iv <- mali_isa$indcls[mali_isa$pval <= 0.05]
pv <- mali_isa$pval[mali_isa$pval <= 0.05]
fr <- apply(mali5_avg > 0, 2, sum)[mali_isa$pval <= 0.05]
fridg <- data.frame(group=gr, indval=iv, pvalue=pv, freq=fr)
fridg <- fridg[order(fridg$group, -fridg$indval),]
print(fridg)## group indval pvalue freq
## QUPR 2 0.8549332 0.003 6
## KAAN 2 0.7534400 0.010 10
## RUSP 2 0.7096551 0.008 7
## SMGL 4 0.6896051 0.012 6# 1 is FallRx
# 2 is MowRx
# 3 is SpringRx
# 4 is Harvest
# 5 is ControlMA_LI MowRx indicators: QUPR, KAAN, RUSP
MA_LI Harvest indicator: SMGL
Finding variation presented by each axis
# axis scores------------------------------
dist1Se <- vegdist(rel_mali5, method = 'bray') # dist in real space
dist2Se <- dist(mali_3d$points[,1], method = "euclidean") # dist. in ordin. space
dist3Se <- dist(mali_3d$points[,2], method = "euclidean")
ax1Se <- lm(dist1Se ~ dist2Se) # distance in real space explained by NMDS axis 1
summary(ax1Se)##
## Call:
## lm(formula = dist1Se ~ dist2Se)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31228 -0.06555 0.01538 0.07752 0.28965
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.61312 0.01304 47.01 <2e-16 ***
## dist2Se 0.17468 0.01567 11.15 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1177 on 251 degrees of freedom
## Multiple R-squared: 0.3312, Adjusted R-squared: 0.3286
## F-statistic: 124.3 on 1 and 251 DF, p-value: < 2.2e-16ax2Se <- lm(dist1Se ~ dist2Se + dist3Se) # distance in real space explained by NMDS axis 1 and 2
summary(ax2Se)##
## Call:
## lm(formula = dist1Se ~ dist2Se + dist3Se)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.238642 -0.053360 0.004834 0.056850 0.241994
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.50292 0.01393 36.11 <2e-16 ***
## dist2Se 0.19466 0.01264 15.40 <2e-16 ***
## dist3Se 0.17845 0.01495 11.94 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09413 on 250 degrees of freedom
## Multiple R-squared: 0.5741, Adjusted R-squared: 0.5707
## F-statistic: 168.5 on 2 and 250 DF, p-value: < 2.2e-16(ax1labSe <- round(summary(ax1Se)$adj.r.squared, 3) * 100) # rsquared of lm of real space explained by axis 1## [1] 32.9(ax2labSe <- round(summary(ax2Se)$adj.r.squared * 100 - ax1labSe, 1)) # r squared of lm of real space explaied by axis 1 and 2 less that of axis 1, to get axis 2 alone## [1] 24.2